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ABSTRACT 

^ ■ Several independent measurements have confirmed the existence of a fluctuation 

excess ((5-Fobs ~ 0.1 nW/m^/sr at 3.6 fim) in the Near InfraRed Background (NIRB) 
up to degree angular scales, whose origin is unknown. By combining high resolution N- 
body/hydrodynamical cosmological simulations and an analytical model, we predict the 
absolute intensity and fluctuations impressed on the NIRB by high-z (z > 5) galaxies 
(some of which harboring Pop III stars, shown to provide a negligible contribution) 
and by simultaneously matching galaxy Luminosity Functions (LFs) and reionization 
j constraints. This strategy also allows us to derive the phenomenological evolution of 

the ionizing photon escape fraction: we find /esc = 1 at z > 11, decreasing to ~ 0.05 
at z = 5. In the wavelength range 1.0 — 4.5 fxm, the predicted cumulative NIRB flux 
is F = 0.2 — 0.04 nW/m^/sr. If galaxies brighter than mi[^ = 28 can be removed, the 
remaining background from stellar sources will still contain a relevant contribution from 
galaxies at z > 5, showing that extraction of the reionization sources signal is feasible. 
However, we find that the radiation from high-z galaxies (including those undetected 
by current surveys) is insufficient to explain such excess: at I = 2000, the fluctuation 
strength is 6F = 0.01 — 0.002 nW/m^/sr, with a relative amphtude 6F/F = 4% almost 
independent of wavelength. Two problems remain unsolved: (a) both the predicted flux 
and fluctuations are considerably lower than the observed excess; (b) the fluctuation 
spectrum is redder (Aq, with p = —1.4) than tentatively measured {p = —3). Both facts 
might indicate that an unknown component /foreground, with a clustering signal very 
similar to that of high-z galaxies, is dominating the NIRB radiation excess we receive 
on Earth. 
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INTRODUCTION 



Observations of high redshift galaxies are essential to understand the reionization of the uni- 

2011al lbh. however, 



verse. The current surveys have reached redshifts ~ 8 — 10 (jBouwens et al. 



201C 



it is generally believed that those detected sour ces, which are usually rare and bright galaxies, are 



not the dominant contributors to reionization ( Choudhury &: Ferraral 120071 : iLorenzoni et alJ 12011 



Jaacks et alJl2012l : lFinkelstein et alJl2012l ). a large portion of galaxies being still below the detection 
limit. Furthermore, it is unlikely to find a substantial fraction of Pop III stars which are believed 
to start reionizatio n in those detected gal axies, since usually they are old and have a long metal 
enrichment history (jSalvaterra et al.ll201lh . 



However, even without detecting galaxies individually, their cumulative flux has the potential 
to provide us key information on the reionization sources. Indeed, the bulk of their radiation 
emitted between the Lyman limit and the visible range, is redshifted into the near infrared band 
at present time. Therefore, the Near InfraRed Background (NIRB) residual, i.e., the remaining 
flux after removal of the contribution from our Solar system, from the Milky Way and the known 
foreground galaxies, would provide us a wealth of information of high redshift galaxies, such as their 
integrated emissivity and their large scale clustering properties, if it is really from those galaxies. 



Measuring the NIRB has a history dating more than two decades (see the review bvlKashlinskv 
2005 ) . Early measureinents gave NIRB >^ 10 nW/mVsr (IDwek fc ArendtJ Il998l : ICorjian et al 



2000; iMatsumoto et al.l l200d : ICambresy et al.l l200ll : iMatsumoto et al.l l2005l ) . It has already been 
confirmed that there indeed exist a non-zero excess when the foreground and the emission from 
the known galaxies in the field of view are removed, and this excess is hard to be attributed to the 



rema ining known galaxy population below the detection limit (jTotani et al.ll200ll : IMatsumoto et al 
2OO5I ). 



However, due to the difficulty in the foreground subtraction (jDwek et al.l l2005l ). it is still 
hardly to accurately determine the absolute intensity, so more recent measurements pay more at- 
tention to intensity fluctuations, which can both exclude the influence of a strong but smooth 
foreground, i.e., the zodiacal light, and help to understand the large scale clustering properties of 



Salvaterra et al 



the s o urces (IKashlinskv et al.ll2002l : iMagliocchetti et al.ll2003l: ICoorav et al 
2OO5I 



2004 



20061) . The recent measurein ents ( Kashlinskv et al. 20051 . 2007a . 



Matsumoto 



et al, 



2012 



Matsumoto et al.l bood . I2OI1I : [Thompson et all l2007al lbl) have obtained the angular power spec- 



tra of the NIRB residual at wavelengths from 1.1 nm to 8 ^m. Those power spectra show that 
the sources of the NIRB residual should have large angular scale clustering properties up to degree 
scales and indicate that the absolute intensity may be far smaller than previous measurements (for 
example, on arcmin scale the measured fluctuations amplitude is about 0.1 nW/m^/sr at wavelength 
3.6 um; for relative fluctuations ~ 10%, the ab solute inte nsity should be on l y ^ 1 nW/m^/sr, see 



e.g. lArendt et al.ll201U : IKashhnsky et al.ll2012l ). Notably, IMatsumoto et al.l (j2005l ) found that the 
intensity and fluctuation excess spectrum are both similar to a stellar spectrum, indicating that 
they might have the same origin. 
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Interpreting the above observations is fairly difficult. ISalvaterra &: Ferraral (|2003l ) and lSantos et al 



(|2002l ) suggested that Pop III stars could possibly be the sources of the residual excess. If true the 
NIRB wo uld be an exqu i site to ol to study Pop III stars and hig h redshift galaxies. However, not 
long after iMadau &: Silkl (|2005l ) and ISalvaterra &: Ferraral (|2006l ) found that this scenario needs a 
very high star formation efficiency, an d may overpredict the high - ^ dropouts galaxiej^. O n the other 
hand, kashlinsky et al] (j2005l . l2007bl ^: iMatsumoto et al.l (l201lh : LKashlinskv et all (|2012l ) still favor 
Pop III stars as the best candidate sources. ICooray et al.l (j2012l ) pointed out that m odels consistent 



with the optical depth of the free electrons after reionization measured by WMAP (jKomatsu et al 



predict much weaker fluctuation ampl i tude than o bserved (by about one order of magnitude) . 
Numerical simulations by iFernandez et al.l ([20101, |2012| ) stressed the importance of the nonlinear 
effects in theoretical calculations to reconcile theory with data. 



This is not th e whole story, however. Some works (iThompson et al.ll2007bl : ICooray et al. 



2007 



Charv et al.ll2008l ) proposed that lower redshift galaxies are possibly the sources of the observed 
fluctuations, but this possibility has become less attractive by now. Indeed, the recent work by 
Helgason et al.l (j2012l ) who reconstructed the emissivity history from the luminosity functions (LFs) 
of observed galaxies, found that fluctuations from the known galaxy population below the detection 
limit are unable to account for the observed clustering signal on sub-degree angular scales. 

At least from a theoretical perspective, one key point is to predict more accurately the NIRB 
contributed by Pop III stars and galaxies before reionization, using models consistent with all 
the current observational constraints, e.g., including high redshift LFs and reionization. If the 
prediction matches the observational levels of both absolute intensity and fluctuations, we can go 
further and investigate in detail the contribution of Pop III stars and high redshift galaxies. If not, 
either our current theories about Pop III stars/high redshift galaxies must be revised, or we may 
need to find alternative explanations for the NIRB residual excess. To construct the star formation 
history from the birth of Pop III stars to the end of reionization as required to reliably predict the 
contribution to the NIRB, we should include the relevant physics of star/galaxy formation, e.g., 
gas dynamics, radiative cooling, supernova explosion and related metal production and transport, 
photoionization and heating, and the complex interactions among these processes. 

In this paper, we attempt to make the most detailed theoretical NIRB modeling so far, enabling 
predictions of both the absolute intensity and angular power spectrum. We calculate the emissivity 
from sim ulations that include all the above physics, and especially a detailed treatment of chemical 
feedback ( Tornatore et al. 2007al lbl). The simulated LFs of high redshift galaxies match remarkably 
well the observations; this is the starting point of our NIRB model. 

The layout of the paper is as follows. In Section [2] we introduce the simulation, and describe 
the steps to calculate the absolute intensity and the angular power spectrum. In Section [3] we 



significant contribution from high-z mini-quasars powered by accretion on to intermediate mass black holes with 
spectra similar to loc al Ultra-Luminous X-r ay sources is disfavored on the basis of the unresolved X-ray background 
intensity constraints ( Salvaterra et al. 20051 ). 
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present our results and compare them with observations. Conclusions are given in Section [H In 
Appendix |A] we compare the different approximate solutions in the analytical c alculation of the 



emissi vity. Throughout this paper, we use the same cosmological parameters as in ISalvaterra et al, 



mm-. Q^=0.26, r?A= 0.74, h=0.73, nb=OMl, n = 1 and (78=0. 8. The transfer function is from 



Eisenstein &: Hul (|l998l ). Magnitudes are in the AB system. 



METHOD 



2.1. The Absolute Intensity 

At z = 0, the absolute intensity of the NIRB observed at frequency z/q is the integrated 
contribution of s ource s whose emission is shifted into a band of central frequency vq. Following 
Salvaterra et al.l (j2006l ). we write it as 

e(i/,z)e-""<=ff(^0'^) 



dz 



cdz 



(1) 



H{z){l + z) ' 

where is the proper distance, is the frequency in the observer frame, = (1 + z)i^o is the rest 
frame frequency, e{i',z) is the comoving specific emissivity, H{z) = ffo\/^m(l + z)^ + is the 
Hubble parameter in a flat ACDM cosmology, c is the speed of light. The effective optical depth 
of absorbers between redshift and z, t^s, is co mposed of two parts: the lin e absorption and the 
continuum absorption; we use the expressions in I Salvaterra fc Ferraral (j2003l ). 



We calculate the emissivity (see also Appendix for further discussion on subtleties related to 
s used in the literature) from the results of the simulations presented 
which include a detailed chemical enrichment treatment developed by 



in 


Salvaterra et al. 


( 


2011 


Tornatore et al. ( 


2007a). 



initial mass function (IMF) ()Salpeterlll955l ). for Pop II stars the mass range is 0.1 — 100 Mq, while 
for Pop III stars the mass range is 100 — 500 Mq. Some recent works indicate that Pop III stars 
may not be as massi ve as several hundred Solar mass as predicted by previous theory, but may be 



limited to < 50 Mq (jHosokawa et al.ll201ll ). Our choice then corresponds to an upper limit to the 



contribution of these sources. Using those simulations, ISalvaterra et al.l (j201ll ) generated the LFs 
of galaxies down to the magnitude far below the current observation limits at high redshifts. In 
the redshift range 5 < z < 10, the simulated LFs match the observed ones almost perfectly in the 
overlapping luminosity range. 

Suppose the specific luminosity of the i-th galaxy in the simulation box is L^,^{z) at redshift z, 
the comoving specific emissivity is 



1 EliLUz) 



e{v,z) 



(2) 



47r V 

where V is the comoving volume of the simulation, N is the total number of galaxies in the 
simulation box at redshift z. 
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In the emissivity calculation, we must correct for rare bright galaxies that a re not caught by 



the si mulation due to the finite box size (10 h ^Mpc). We follow the steps in ISalvaterra et al 



(|201ll ). We first calculate the absolute magnitude corresponding to the mean luminosity of the two 
brightest galaxies in the simulation box, Muv,up- The contribution (to be added to the numerator 
in Eq. ([2|)) from galaxies brighter than this magnitude is obtained by integration 

A'/uv.up 



Lr'iz) = V 



Ll{z) 



25 



$(Muv,2)dMuv, 



(3) 



where L\^{z) is the luminosity of the brightest galaxy in the simulation (we assume all rare bright 
galaxies have the same Spectral Energy Distribution (SED) of this one), Luv is the luminosity cor- 
responding to the UV absolute magnitude Muv- The wavelength used to calculate the absolute UV 
magnitude in this paper is 1700 A, z^uv is the frequency corresponding to this wavelength. In ob- 
servations, the selected wavelength corresponding to the UV abs olute magnitude may be somewhat 



different in different meas urements and at different redshifts (jBouwens et al.l 120071 : lOesch et al 
20ld ; iBouwens et al.ll2010l ). however, our results are not sensitive to such di fferences. For t he LF 



$(MuV)'Z) in the redshift range 5 < z < 1 0, we use the Schechter formula (jSchechterlll976l ) with 
the redshift-dependent parameters given bv lBouwens et all ()2011bl ) (see their Sec. 7.5), who fitted 
the observed LFs in z ~ 4 — 8 and extrapolate them to higher redshifts. For redshifts above 10, 
we simply add an exponential tail normalized to the simulated LF amplitude at Muv, up- We find 
that the entire correction has only a small impact on our results. As discussed below (see also right 
panel of Figure [3]) , ~ 90% of the high-z galaxy contribution to the NIRB intensity comes from 
sources at 5 < z < 8 where the correction is at most 12%. This correction is also applied to the 
calculation of ionizing photons below. 

For each galaxy, the radiation comes from two different mechanisms: the stellar emission 
and the nebular emission. The former comes directly from the surface of stars, while the latter 
is generated by the ionized nebula around stars and depends on the fraction of ionizing photons 
that cannot escape into the intergalactic medium (IGM), i.e., 1 — /eso where /esc is the escape 
fraction. Ionizing photons escaping from galaxies would ionize the IGM; such ionized gas could 
also pro duce the nebular emiss ion. However, due to th e very low recombination rate, as shown 
in, 



e.g., 



Nakamoto et al.l (|200ll ) and ICoorav et al.l (|2012l ) , its emissivity is much weaker than the 



radiation from galaxies, so we ignore this contribution in th is paper. The IGM contribution to the 
NIRB fluctuations is also negligible (jFernandez et al.ll20ld ). 



To determine the escape fraction averaged over the galaxy populations present at a given 
redshift, we proceed as follows. First we compute the number of ionizing photons emitted per 
baryon in collapsed objects as: 



1 



E 



N 



N 

E 



1=1 ^'^gas i=l 



L^IIi'^ILi r7i\ ii'rll i H i I ^III 11 H 



(4) 



where is the emission rate of ionizing photos from Pop II stars (this quantity depends on both 
the age and the metallicity of the stellar population) corresponding to a continuous star formation 
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rate 1 Mfr^yr 



We derive this quantity from the Starburst99 template^ ( Leitherer et al. 1999 : 
Vazquez Sz Leithereiil2005l : iLeitherer et al.ll20ld ) adopting the mean age, r^^'*, and metalhcity, Z^, 



of each simulated galaxy, is the emission rate of ionizing photons for Pop III stars according 



to 



Schaererl ^200^ ) @. Mg^g is the gas content, ' is the mean star formation rate of Pop II 



stars in this galaxy, while M^^^'* is the cumu lative mass of Pop III stars. We use a mean lifetime 



= 2.5 X 10^ yr for massive Pop III stars ()Schaererll2002l : ISalvaterra et al.l 



2011 



We then compare the above quantity with the number of ionizin g photons per bary on in 
collapsed objects, A'^ion, required by interpreting the observations as in iMitra et all (j2012al ) (the 
"mean" value) to get the escape fraction, i.e., /esc = min( '^^°" ,1.0), where C is the clumping 
factor. Throughout this paper we assume C = 1 to get the minimum /esc therefore the maximum 
contribution of the nebular ei nission to the NIRB. We note that th e clumping fac tor could be highe r 



than 1 even at high redshifts (jPawlik et al.ll2009l : iShull et al.ll2012l ). For example. IShull et all (|2012l ) 
gives C ~ 3 (1.7) at z = 5 (9) by numerical simulations. Our nebular emission is therefore reduced 
by about 10% to 60% from z = 5 to z = 9 if this clumping factor is adopted. However, for Pop II 
stars which are the dominant contributors to the NIRB, the nebular emission is much smaller than 
the stellar emission. So the final reduction in the NIRB would be much smaller. Furthermore, we 
will show later that the flux from high redshift galaxies and Pop III stars is not able to explain the 
observed fluctuations level, so that a reduction in the NIRB intensity would in any case strengthen 
this conclusion. We plot the derived /esc as a function of redshift in Figure [H There is a clear trend 
of an increasing escape fraction towards higher redshifts; it reaches 1 at z ~ 11. At z = 5, the final 
redshift of the simulation, /esc ~ 0.05. Although required by reionization data, an increasing trend 
of /esc (-2) has not yet been fully understood theoretically in spite of the several, often conflicting, 
studies on this problem. 



Based on observations. Ilnoue et al.l (120061 ) concluded that /esc > 0.1 when z > 4; by combining 
the observations of Lyman a absorption and UV LF, and also using N-body simulation s and semi- 
analytical prescriptions to model the ionizing background. ISrbinovsky Sz Wyithd (120101 ) found that 
for galaxies at z ~ 5.5 — 6, if the minimum mass of star f ormin g galaxies corresponds to the 
hydrogen cooling threshold, /esc ~ 0.05 — 0.1; IWvithe et al.l (j2010l ) used the star formation rate 
derived froin gamma-ray burst observations to conclude that in the redshift range 4—8.5, /esc ^ 0.05; 
Wise (fe CenI (j2009l ). by radiation hydrodynamical simulations, found that at redshift 8 for galaxies 
with Mvir < 10^-5 Mq, /esc ~ 0.05 - 0.1, while for more massive g a laxies /esc ~ 0.4, if a normal 
IMF is adopted; also via simulations, iRazoumov &: Sommer-LarsenI (120101 ) found /esc ~ 0.8 when 
z = 10. The escape fraction derived by us is broadly consistent with those values. The important 
difference however is that our derivation of /esc matches both the LF and the reioniz ation history 



simultaneously, i.e., a more phenomenological derivation. Different from our approach. iMitra et al 



^http://www.stsci.edu/science/starburst99/docs/default.htm 

■^Note the different dimensions of and q™: the former corresponds to per unit star formation rate while the 



latter corresponds to per unit stellar mass. 



-7- 



(|2012bl ) computed the LFs of high redshift galaxies by means of semi-analytical models and derived 
the star formation efficiency required to match the observed ones. They found /esc ~ 0.07 at 
z = Q and /esc ~ 0.16 at z = 7, which are consistent with our /esc ~ 0.06 (0.18) at those two 
redshifts. At higher redshifts, however, their escape fraction is somewhat lower than ours. 

As a final remark, we underline that when computing the escape fraction, we do not make 
a distinction between Pop III and Pop II stars. Since in the calculation of f-^N^ we account for 
ionizing photons from both of them, one can think of /esc as a kind of "effective" escape fraction 
averaged over the galaxy population. In principle, /esc for Pop III stars should be higher due to their 
harder spectrum. However, as we will see in Section [3l Pop III stars only contribute a negligible 
flux to the present-day NIRB; a more detailed modeling is then not necessary. 




10 
z 



15 



Fig. 1. — Escape fraction evolution from joint LF-reionization constraints. 



The luminosity of the i-th galaxy, L]^, is the sum of the contribution of Pop II and Pop III 
stars. For Pop II stars we use the age and metallicity dependent spectrum templates provided 
by the Starburst99 code. The nebular emission contribution has been renormalized by adopting 
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the escape fraction computed above. In addition to free-free, free-bound and two-photon emission 
which are already included in Starburst99, we add the Lyman a emission to the template by using 
([Fernandez fc Komatsull2006l ) 

l^{y,T^^\Z\z) = fahpu^c^iu - z.„)(?g(r"^Z^)[l - f,,,{z)], (5) 

in which = 0.64 ( Fernandez fc Komatsul |2006| ) . hp is the Plank constant, z/q = 2. 47 x 10^^ Hz 
i s the frequency of Lyman a photons. We use the line profile (^(z^ — z/q) provided in I Santos et al 

(|2nn2l ): 







l'a)'^exp[-Ui,{z)/\u - Ua 



where 



u^{z) = 1.5 X 10 



11 



0.019 



h 

07 



[1 + 



when f < fa 

when U > Ua: 

Hz 



is the fitted form of results given in lLoeb &: Rvbickil (jl999l ). 



(6) 



(7) 



For the template of Pop HI stars, l]}^, we still use the spectrum in ISchaereil (|2002l ). but 



renor malize the nebular em ission part by the factor 1 — /esc- The luminosity of the i-th galaxy is 



then (Salvaterra et al. 



2011 



the Lyman a emission in Eq. ([5]) is already included in l]} here. 



,111/ 



rIII,j^III 



(8) 



Given the above luminosity for each galaxy, we then get the emissivity according to Eq. ([2]). 
As an example, we plot the ve{v,z) at redshift 12.0, 9.0 and 6.0 respectively in Figure [2j At high 
redshifts, the escape fraction is ~ 1.0, yielding a very weak Lya line, since such emission is produced 
by recombinations of the ionized nebula around stars, which is ionized by the fraction (1 — /esc) of 
the ionizing photons (dash-dotted line in Figure [2|). At lower redshifts, due to a decreased escape 
fraction, the Lya emission is clearly seen in the spectrum, as shown by the dashed and solid lines. 

By focussing on the region of most interest here (energies below 10.2 eV), we find that spectrum 
becomes increasingly flatter with time. For example, at z = 12, the slope of z^e(i^, z) oc with 
/3 ~ 2, while atz = 5/3~1.2. This is clearly the result of an aging effect enhancing the rest frame 
optical/IR bands flux with respect to UV ones. Since the NIRB from z > 5 galaxies is dominated 
by lower redshift galaxies (5 < z < 8), we do not expect to get very steep NIRB spectra, as indeed 
we show in Sec. [3] 



2.2. NIRB Fluctuations 



The angular powe r spectrum of the fluctuations of the specific intensity field given by Eq. ([T]) 
is (jCoorav et al.lboid l 

dz dr 



r2(z)(l + z)2dz 



=(z.,z)e--^«('^0'^)]2Pgg(A:,z), 



(9) 
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Fig. 2. — The emissivity of simulated galaxies at redshift 12.0 (dash-dotted), 9.0 (dashed) and 6.0 
(solid) respectively. 



where the comoving distance is r{z) = ' ^ ~ However, in ob servations what is 



really measured is the fluctuations of the "flux", which is defined as F = vqI^^ ([Kashlinskv 



2005; 



Thompson et al.ll2007al ). To be consistent with this convention we replace e{v,z) with VQe{v^z) 



iye{i',z)/{l + z), the above formula is then reduced to 



C, 



cdz 



[z^e(i/,z)e-^-ft('^"'^)]2 



Pgg{k, z); 



(10) 



i/(z)r2(z)(l + z)4 

without danger of confusion we ignore the superscript. Eq. ([TO]) assumes that the luminous proper- 
ties of galaxies are independent of their locatio n, so the only fact or that determines the fluctuations 
of the NIRB is their spatial fluctuations (see lShang et al.ll2012l for an improved model), which is 
described by the galaxy-galaxy power spectrum, Pgg{k,z). 



The 10 h-^ Mpc box size of ISalvaterra et al.l (|201ll ) simulations is too small to provide us 
with the large scale correlation function of galaxies (for sources at z = 6, the comoving trans- 
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verse separation corresponding to 1° angular size is about 100 /j-^Mpc), so we use the halo model 
(jCoorav &: Shethll2002l : ICooravll2004l ) to calculate the galaxy-galaxy power spectrum. This power 
spectrum is composed of two parts, the one-halo term which means contributions from the corre- 
lation of galaxies (including the central galaxies and satellite galaxies) in the same halo, and the 
two-halo term from galaxies in different halos: 



Pgg{k, z) — i^gg (/c, z) + -Pgg (/c, z) 



(11) 



where, by assuming that the distribution of galaxies in a halo traces the profile of dark matter, and 
the mean number of central galaxies and satellite galaxies in a halo with mass M are {Nsax) and 
(A^cen), respectively. 



"'"'"^'^dM— ^^^^^^^^ 
dM 



n 



gal 



and 



P^^{k,z) = PUk,z) 



dM b{M, z)}-^E±Jl±^ 

,{z) dM ngai 



-u{M, k) 



(12) 



(13) 



In above expressions, M^i^{z) is the minimum mass of halos that could host galaxies, i.e., the 
minimum mass of halos that contribute to the emissivity. We set it to be the minimum mass of 
halos that contain stars in our simulations, which is ~ (2 — 8) x 10 Mq, depending on redshift. 
Mraax{z) is the maximum mass contributing to the emissivi ty and the clustering ( we describe 



how we determine it later on); dn/dM is the mass function (jSheth &: Tormen 



199S 



Sheth et al 



200lh. while u(M,k ) is the normalized Fourier transform of the halo profi l e. Fo r a NFW profile 
(jNavarro et al.lll997l ). the analytical expression is given bv lCoorav &: ShethI (j2002l ). and we use the 
concentration parameter, cm, from iPrada et al.l (j2012l ) which fits simulations well. However, we 



find th at our results are insensitive to the use of different concentrations, as e.g., from lZehavi et al 



(j201ll ). Even if a concentration largely different from the above ones is adopted, its effects are 
non- negligible only in the one-halo term which dominates the signal at small scales, where galaxy 
clustering is well below the shot noise. Therefore, since we are mainly interested in the large-scale 
(> 1') clustering, our conclusions are unaffected by the adopted value of cm- Finally , 5(M , z) is the 
bias of halos relative to the dark matter; we use the formula given by lTinker et al.l (l20ld ) and the 
fitted paraineters f rom simulations therein. This formula predicts higher bias for massive halos than 
Sheth et al.l (120011). but i s mor e close to simulations. The linear matter power spectrum, Pi[^{k, z), 
is from lEisenstein &: Hul (|l998l ). 



The mean number of central galaxies a nd satellites in a h alo with mass M is modeled by the 
halo occupation distribution (HOD) model (jZheng et al.ll2005l ). 



(^cen) = \ 



1 + erf 



logioM - logioMn 



and 



sat/ 



1 + eif 



0"logioJ\/ 

logipM - logio2Mmii 
O'logloA/ 



M 
M^i 



(14) 
(15) 
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We adopt the parameters Mgat = IS Mmim o'\r,^,nM = 0- 2 and as = 1.0, w hich are from bot h 



simulations and semi-analytical models (jZheng et al.ll2005l ). and observations (jZehavi et al.l 12011 



With the mean number of central and satellite galaxies in each halo, the galaxy number density is 
simply 

/ dM—{{N,,,,) + {Nsat)). 



n 



gal 



dM 



(16) 



In addition to the above galaxy clustering, Poisson fluctuations in the number of galaxies 
would generate shot noise in observations, whose power spectrum dominates at small scales. If the 
redshift derivative of the number of sources with flux between 5 and S + dS is the angular 

power spectrum of such shot noise is 



Cl 



SN 



1 

AO 



dz / dSS^ 



d^N 



1 



dSdz An 



dz / dMS'^ 



d^N 
dMdz' 



where Afi is the beam angle. Considering that 

d^N dn 



, dr 



and 



dMdz dM "' dz' 

O - 



47rr2(l+z) ' 

where Ly{M) is the luminosity of halos with mass M, the shot noise power spectrum is 



dz 



g-TcffCfcz) dr 
r2(l + z)2dl 



dn 
dM^ 



-dM. 



(17) 

(18) 
(19) 

(20) 



As we assume the luminous properties of galaxies are independent of their location, we can simply 
use an average light-to-mass ratio that is independent of halo mass. 



L^{M) 
inM 



Ly{M) 



dn 
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(21) 



As in the galaxy clustering, by replacing e(z^, z) with z^e(i/, z)/(l + z), we finally get the shot noise 
angular power spectrum 



C, 



SN 



cdz 



H{zy{z)il + zY- 



where 



i'e{v, z)e '^^ff 



Ph 



dMM^ 



dn 
dM' 



(22) 
(23) 



dMM^. 
dM 



and the halo mass density ph 

In observations, galaxies are generally removed down to a certain limiting magnitude, mum, to 
get the NIRB residual fluctuations. Thus, we need also to remove galaxies brighter than this limit 
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in the simulation box and in the bright-end. At A n = 1-6 //m iThompson et alj (j2007al ) adopted 
m]jm — 28; at An = 2.4 fj,m iMatsumoto et all (j201ll ) reached ~ 23; and at An = 3.6, 4.5 fim 
Kashhnskv et al.l (|2012l ) used mnm — 25. The apparent Hmiting magnitude (at wavelength An) is 
converted into the rest frame absolute magnitude (at wavelength Ao/(l + z)) Mx„/(^i+z) by 



M 



\o/{i+z) = m-im - DM{z) + 2.51ogio(l + z) 



(24) 



where DM{z) is the distance modulus (jHelgason et al.ll2012l ). By using a light-to-mass ratio con- 
structed from the simulation, we determine the maximum halo mass M^a.-K{z) in the calculation of 
galaxy-galaxy power spectrum according to the above absolute magnitude. Things are slightly more 
complicated when calculating the absolute intensity and the spectrum of fluctuations, since they 
both depend on wavelength , while the limiting magnitude in observations at different wavelength 
is different. In this case we simply give the theoretical prediction without removing any sources in 
the simulation box and the bright-end, as shown by Figure [3] and Figure [71 We use Mmax = oo and 
then we discuss the effects of galaxy removal, i.e., in Figure HI Throughout this paper we adopt 
Zmin = 5 and Zmax = 19 uuless otherwise specified. 



3. RESULTS 

We start by presenting the absolute intensity of the NIRB observed at ^; = in the (observer 
frame) wavelength range 0.3 — 10 fim using Eq. ([2]) and the bright-end correction as defined by Eq. 
([3]). Figure [3] (left panel) shows the predicted intensity i^qIuq (in cases without risk of confusion, we 
will refer to this variable as "flux" ) when all sources with z > 5 are included; also shown separately 
are the contributions from Pop II and Pop III stars. The flux peak value is 0.2 nW/m^/sr at 
Aq = 0.9 fim, and decreases to 0.04 nW/m^/sr at An = 4.5 ^m. The small bump on the left of the 
peak is due to intergalactic Lya absorption by intervening neutral hydrogen. 

We find that in our case, the Pop III contribution is almost negligible (it never exceeds 1%). 
This results from the fact that in the simulation the Pop III star formation rate is about three 
orders of magnitude lower than that of Pop II stars at z = 10; the ratio is even smaller below this 



redshift (jTornatore et al.ll2007bl ). Stated differently, halos with the highest Pop III stellar fraction 



are usually smaller ari d less luminous, and their contribution to the total luminosity is very low 



(jSalvaterra et al.ll201ll ). So one should not be surprised that Pop III stars contribute negligibly to 
the NIRB. This prediction makes it more difficult to find Pop III signatures by means of NIRB 
observations. 

In the right panel of Figure O we plot the contribution from sources above redshift 5.0, 8.0 
and 12.0, respectively. It is clear that contributions from sources located at 5 < z < 8 dominate, 
providing about 90% of the flux from all sources with z > 5. Most of those sources are the low- 
luminosity galaxies that cannot be detected individually in current surveys: they are believed to 
be the major contributors to reionization. In principle, then, the NIRB could be a perfect tool to 
study the reionization sources without detecting them individually. 
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Fig. 3. — NIRB flux from high-z galaxies, vqIuq, at 2; = as a function of wavelength in the observer 
frame. Left panel: Contributions from Pop II (dashed line), Pop III (dash-dotted) stars, and their 
sum. Since the contribution from Pop III stars is very small, the solid line and the dashed line are 
almost identical. Right: Contributions from sources in different redshift ranges: z > 5 (solid line), 
z > 8 (dashed) and z > 12 (dash-dotted). In b oth panels, circles (squar es) refer to the measured 
isotropic emission by iMatsumoto et al.l (j2005l ) (jWright &: JohnsonI I2OOII ) after subtraction of the 
contribution from galaxies with z <b. 



Using the LF constructed in lHelgason et al.l (j2012l ). we calculate the integrated light of galaxies 
with z < 5. We then plot in Fig ure [3l the remaining flux after subtrac ting this component frorn 



the isotropic emission measured in IMatsumoto et al.l (j2005l ) (circles) and IWright &: JohnsonI (j2001 



(squares). Figure [3] clearly shows that even after such subtraction, the remaining flux in the 
measurements is > 100 x larger than the theoretical prediction, indicating the need for another yet 
unknown component. 

We plot the NIRB flux before and after removal of bright galaxie s in Figure HI at wavel engths 
from 1.25 /im to 4.5 /im, corresponding to the band from J to M in iHelgason et al.l (j2012l ). The 
crosses refer to flux from all galaxies with z > 5 in our work, while t he solid line correspond to 
the flux from all galaxies in the default model of IHelgason et al.l (j2012l ). After removal of galaxies 
brighter than nin^ = 28, the fl ux from z > 5 galaxi es in our work is shown by diamonds, while 
flux from remaining galaxies in IHelgason et al.l (j2012l ) is shown by the dashed line. The flux from 
all galaxies (solid line) is about 1-2 orders of magnitude larger than that from galaxies with z > 5 
(crosses) in our work. Hence, without bright galaxy removal, z < 5 galaxies largely dominate the 
NIRB flux. However, if we remove the galaxies down to = 28, the flux from the remaining 
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galaxies at all redshifts in iHelgason et al.l (|2012l ) (dashed line) is comparable to that from the 
remaining galaxies with z > 5 in our work. Even considering the uncertainties on the faint-end of 
LFs, in the source-subtracted flux, galaxies at z > 5 still contribute at least ~20%— 30%. So we 
find that at least in principle we can indeed access the signal of reionization sources by subtracting 
the bright galaxies. 

The fluctuations of the NIRB after subtracting galaxies down to the detection limits of obser- 
vations, + l)Q/(27r), at Aq = 1-6, 2.4, 3.6, 4.5 //m are shown by the thick solid line in each 
panel of Figure [Sj the shot noise is also shown by a thick dashed line in each panel. At large scales, 
i.e., I < 10^, the galaxy clustering is larger than the shot noise level, implying that z > 5 galaxies 
can indeed produce large scale clustering signatures. So the NIRB fluctuations do have the poten- 
tial to provide us wealth of information of the clustering of the undetected reionization sources. 
However, compared with the o bserved values (shown by filled circles with errorbars) which are 



usually at 0.1 nW/m^/sr level ( Thompson et al.ll2007al : iMatsumoto et al.ll201ll : iKashlinskv et al, 
20121 ). our predictions are only ~ (2 — 4) x 10"'^ nW/m^/sr, i.e., much smaller than observations, in- 



dicating again the existence of unknown component (s) we are missing. Somewhat surprisingly (but 
interestingly) the missing component has a clustering signal very similar to that of high redshift 
galaxies and extending to degree angular scales. Obviously, this component /foreground should be 
removed before we are ready to exploit the NIRB to study reionization sources. 

In each panel we also plot the galaxy clustering and the shot noise in the case without removal of 
galaxies by thin lines. Except for the 1.6 /im case in which galaxies are removed down to mnm = 28, 
at different wavelengths the large scale galaxy clustering is almost unchanged, no matter whether 
the bright galaxies are removed or not. This implies that the corresponding limiting magnitude is 
not deep enough to see obvious changes of the galaxy clustering. 



Next, we define 6F = + l)C//(27r) as the fluctuation amplitude and plot its ratio to F = 
i^oIuQ in Figure [6l Such relativ e fluctuation is almost independent of the observed wavelength, Aq 
(see also [Fernandez et al.ll20ld ). and 6F/F ~ 4% at / = 2000, with the only exception of the 1.6 /im 
band where it is somewhat lower than in redder bands, as a consequence of a deeper (mwrr , = 28 ) 



gala xy removal. N i cely, the relative fluctuation agrees with that found bv [Fernandez et al.l (j2010l ) 
and ICoorav et al.l ()2012l ). In addition, 5F/F increases with Zmin, i-e., high redshift sources have 
higher relative fluctuations. For example, 5F/F = 7% for Zmin = 8.0, while it reaches 12% if 
Zmin = 12.0 is adopted: although both 6F and F decrease for larger Zmm, the P^^{k, z) term in 
6F reduces the decline rate of 6F, since k = l/r{z) corresponds to larger scales at higher redshifts 
if / is flxed. In other words, the higher is the redshift of the sources, the larger is the scale (and 
consequently P'^g) contributing to the clustering signal. The relative fluctuat ion SF/F is only 



weakly dependent on the intrinsic properties of galaxies (jFernandez et al.ll2010l ). but more so on 
the spatial clustering feature^. Thus, 5F/F is a key indicator to identify NIRB sources; yet, in 



^As a sanity check, models that find a fluctuation level matching observations should avoid over-predicting the 
absolute intensity at the same time. 
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practice, it is hard to get an accurate absolute intensity. 

The spectrum of the fluctuations, 6F{Xq) from all galaxies with z > 5 at I = 2000, shown 
in Figure [71 has a slope Aq, with p = —1.4 above 1 /um. Such slope is essentially the same 
as that of the flux, reflecting the above mentioned wavelength independency of 5F/F. There is 
circumstantial evidence from observations that the spectrum is instead "bluer" '^g^' rern iniscent 



of the Rayleigh- Jeans regime of a blackbody spectrum, e.g. iMatsumoto et al.l (|2005l . l201ll )). or at 
least steeper than —1.4. Our spectrum is redder basically for three reasons. First, the spectrum of a 
galaxy is the combination of spectra of stars spanning a wide mass range. Usually, it is redder than 
a typical massive star, unless a very top heavy IMF form is adopted. Second, galaxies close to the 
end of reionization {z ^ 5 — 8) contribute most to the integrated NIRB from high- 2; galaxies. Due 
to their larger ages and higher metallicity, those galaxies are usually redder than earlier systems. 
Third, Starburst99 includes not only the main sequence but also later evolutionary stages, during 
which stellar spectra become redder. The wavelength range 1.0—4.5 /im corresponds to a rest-frame 
energy 2 — 9 eV provided the redshift of sources is ~ 6. The slope of the fluctuation spectrum we 
predict is consistent with that of the flux spectrum around this redshift, i.e., the combination of 
the speciflc luminosity of galaxies with different ages, metallicity and so on in this energy range. If 
the data are confirmed, an important tension remains, as the blue colors in observations can hardly 
be attributed to the galaxy population consistent with the observed LFs and used in this paper to 
predict the NIRB. 



4. CONCLUSIONS 

By combining high resolution N-body/hydrodynamical cosmological simulations and an ana- 
lytical model, we have predicted the absolute intensity and fluctuations impressed on the NIRB by 
high-z (z > 5) galaxies, some of which harboring Pop III stars. This is the most robust and detailed 
theoretical calculation done so far as we simultaneously match also the LFs and reionization con- 
straints. The simulations include the relevant physics of galaxy formation and a novel treatment 
of chemical feedback, by following the metallicity evolution and implementing the physics of Pop 
Ill/Pop II transition based on a critical metallicity criterion. It is able to reproduce the observed 
UV LFs over the redshift range 5 < z < 10, and extend it to faint magnitudes far below the 
detection limit of current observations. 

We directly calculate the stellar emissivity from the simulations. We use Starburst99 to gen- 
erate metallicity and age dependent SED templates, then calculate the luminosity for each galaxy 
according to its current star formation rate, stellar age and metallicity, instead of using a constant 
metallicity and average main sequence spectrum template. Except for the mass range of the IMF 
that has already been fixed in the simulation, there are no other free parameters in the calculation 
of the emissivity. 

By comparing the number of ionizing photons produced per baryon in collapsed objects, fi,N-y, 
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in the simulation and the ionizing photon rate A'ion ~ /esc/* -^7 deduced from data-constrained 
reionization models, we get the evolution of the escape fraction of ionizing photons, /esc(-z)- We 
find /esc ~ 1 at z > 11, decreasing to ~ 0.05 at z = 5. This escape fraction is used to renormalize 
the nebular emission of Pop III and Pop II stars in the emissivity. 

Pop III stars are unlikely to be responsible for the observed NIRB residual, as their contribution 
is very small (< 1%) of the absolute intensity in our calculation. This is the natural result of the 
much lower star formation rate of Pop III stars compared with Pop II stars in the simulation, since 
even metals from a single Pop III star c ould enrich above the critical metallicity a large amount 



of gas around it (jTornatore et al.ll2007bl ). The formation of Pop III stars is regulated by such a 
chemical feedback mechanism, which limits their contribution to the NIRB. However, a rapid Pop 
Ill-Pop II transition brings al so a little advantag e in terms of integrated emissivity, due to the 



longer lifetime of Pop II stars (ICooray et al.ll2012l ). 



We predict that in the wavelength range 1.0 — 4.5 ^m, the NIRB flux F = vqIv^ ~ 0.2 — 
0.04 nW/m^/sr, while the fluctuation strength is about 5F = 0.01 — 0.002 nW/m^/sr at / = 
2000, when all z > 5 galaxies (and their Pop III stars) are accounted for. If we remove galaxies 
down to m\\m = 28, the above flux level is only slightly reduced; however, by comparing with 



Helgason et al.l (|2012l ). we find that the flux from z < 5 dramatically decreases and the remaining 
becomes comparable to the predicted signal. This implies that in principle it is possible to get the 
signal from reionization sources by subtracting galaxies down to a certain magnitude. 

The relative fluctuation amplitude, 5F/F, at / = 2000 is ~ 4%, almost independently of the 
wavelength. This ratio may be helpful to investigate the clustering features of the NIRB sources, 
since the intrinsic properties of galaxies almost cancel out. Despite the difficulties in measuring 
the absolute intensity accurately, it could be treated as a quality indicator in the data reduction 
process: if a much higher/lower ratio is obtained from the data, this might suggest that a more 
careful analysis work is required to extract the genuine contribution from reionization sources. 

In spite of being accurate and consistent with LFs and reionization data, thus offering a 
robust prediction of the NIRB contribution of high-z galaxies that likely reionized the universe, our 
interpretation leaves two open and puzzling questions. 

First, both the predicted absolute flux and the fluctuations are considerably lower than ob- 
servations, indicating that in addition to the contribution from the expected high-z galaxy pop- 
ulation (and Pop III stars), we should invoke some other - yet unknown - missing component (s) 
or foreground (s) dominating the currently observed NIRB residual. The missing component has a 
clustering signal very similar to that of high redshift galaxies and extends to degree angular scales. 
Obviously, this component/foreground should be removed before we are ready to exploit the NIRB 
to study reionization sources. On the other hand, sources located at 5 < z < 8 provide about 90% 
of the flux from all sources with z > 5 in our simulation; these are the faint galaxies currently 
undetected by deep surveys. Thus, if the above mentioned additional spurious sources/foregrounds 
can be removed reliably, the NIRB will become the primary tool to investigate the properties of 
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the reionizing sources. 



Second, the fluctuation spectrum is redder /flatter (Aq, with p = —1.4) than that suggested 
by the data, which seems to favor a slope close to a black body spectrum in the Rayleigh-Jeans 
regime (p = —3). In principle this could only be explained in the rather unlikely case in which the 
spect rum is dominated by massive Pop III stars with a very peaked IMfi. Our study and related 
ones (|Maio et al.1120101 : IWise et al.ll2012l ) however strongly disfavor such an extreme scenario. As 
"normal" galaxies observed in the HUDF tend to have SEDs consistent with our predictions, the 
blue color of the spectrum, if confirmed, might again point towards the fact that an unknown 
component /foreground is dominating the residual NIRB radiation we receive on Earth. 
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AN ANALYTICAL DERIVATION OF THE EMISSIVITY 



At redshift the emissivity (defined as energy emitted per unit comoving volume per unit 
frequenc y per unit steradian here) of stellar population at frequency v is given by the following 
integral (jFernandez &: Komatsul 120061 ): 



e(i/,z) 



Y rm-2 



Air /., 



Lu{rri)ni,[rri)dm, 



(Al) 



■mi 



where Lu{m) is the specific luminosity of a star with mass m, mi is the minimum mass of stars 
while 1712 is the maximum mass, ni,(m) is the number density of alive stars between m and m + dm, 
which is written as 

l-t(z) 

n^{m)= I h^{m,t')dt', (A2) 

Jt{z)-T{m) 

where t{z) is the age of the universe at redshift z, T{m) is the lifetime of a star with mass m. For 
Pop II stars with metallicity 1/50 Zp^, us e ful fit ting formulae for these quantities as a function of 
m are collected in [Fernandez &: Komatsul (j2006l ). Eq. ()A2p means that only stars formed between 
t{z) — T{m) and t{z) can still emit photons at time t{z). The formation rate of stars with mass 
between m and m + dm, h^{m, t'), is 



n,(m,t') = ^f{m), 



(A3) 



mf{m)dm, 



(A4) 



/■m2 

in which f{m) is the normalized stellar IMF, i.e., / f{m)dm = 1 and 

J mi 

while 

PV(0 = A^ / M-^iM,t')dM 
Jm,,,;,, dMdt' 

is the comoving star formation rate density in halos above mass Mmim provided a fraction of 
baryons are converted into stars. 

Two approximate solutions can be found under particular circumstances. If the star formation 
rate is almost constant over the time interval T(m), i.e., r(m) < t^Y-{z), where the star formation 
\ Pi,iz) 

time scale tsFi-^j ~ 



then 



Jt(z)- 



t{z)—T{m) 



p^{t')dt' « p^[t{z)]T{m), 



and the emissivity is approximated as (jFernandez &: Komatsu 



2006 



Fernandez et al.ll2010l ) 



e(i/, z) = — ^ [ L,y{m)T{m)f{m)dm, 



Air m^ 



(A5) 



mi 



which is usually used for relative massive stars with short lifetime, we name it 'Approximation 1" . 
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On the other hand, if r(m) is longer than the age of the universe (this is true for stars of 
smaher mass, and means that no stars die). 



dt' 



t{z)-T(m) 



t{z) 



dt' 



(A6) 



the emissivity becomes ([Fernandez &: Komatsul 120061 ) 

1 P.{z) 1"^^ 



€{U,Z 



Ly{m) f {rn)dm. 



(A7) 



mi 



This also holds true if T(m) is much longer than the star formation time scale tsrC^), i.e. the death 
of stars is less significant compared with the formation of new stars, so that 



t(z) 



dt' 



(A8) 



lt{z)-T{m) 

and the emissivity could also be approximated as Eq. (IA7p . and we call it "Approximation 2". 

In reality, a galaxy is composed of stars with different mass; some of them may have lifetime 
longer than tsF, while others not. In this case a "Hybrid" approximation could be used. 



'p*{z) 



mt 



Ly{rri) f [vn^jdm + 



P*iz) 



mi 



m2 



L^{m)T{m)f(m)dm 



nit 



(A9) 



where mt is the stellar mass determined by the condition Turrit) = tsF- 

We compare the full analytical solution of Eqs. ()AlllA4p with these three approximate so- 
lutions. Furthermore, we will also consider the emissivity obtaine d by adopting the Starburst9 9 
template at Z = 1/50 Zq instead of the simplified fitting formula bv [Fernandez &: Komatsul (j2006l ) . 
In this case the emissivity is given by 



e{i',z) 



1 

47r 



tiz) 



Lu,SB99[z,t{z) - t']p^{t')dt' , 



(AlO) 



where -L^,sB99 is the luminosity per unit mass (note that here for integration purposes we use the 
burst star formation model) from Starburst99. 

In our work, the mas s range of Pop II stars is 0.1 — 100 Mp), while t he fit ted f ormula of the 



main sequence age used in [Fernandez &: Komatsul (|2006l ) . iFernandez et al.l (|2010l ) and lCoorav et al 



(|2012l ) (taken from ISchaered |2002| ) is based on data of massive stars. To avoid introducing more 
uncertainties, in this comparison we adopt a mass range 1 — 100 Mq for Pop II star s. We checked 



that f or Pop II stars with mass 1 Mq, the fitted main sequence age still agrees with lGirardi et al 
(|200d ). 

Since Pop II stars are found to contribute much more than Pop III stars to the NIRB (see 
Figure [3]), and stellar emission is the dominant component, we neglect here the nebular emission. 
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In this case L,^ is .just a bla c kbody spectrum, and we truncate it at hv = 13.6 eV, as done in 
Fernandez fc Komatsul (|20od '): IPernandez et alJ (|20ld ^: ICoorav et ali ^20li ). 



We plot the emissivity at z = 10 calculated by different methods in Figure[8]for a star formation 
efficiency f^, = 0.01 and a minimum mass Myain = 10^ Mq. It is not surprising that "Approximation 
2" overestimates the emissivity, since it assumes no stars die. However, we find that for the stellar 
mass range 1 — 100 Mq, "Approximation 1" also overestimates the emissivity when hv <S eV. This 
happens because of the contribution of low mass stars whose lifetime is even longer than the age 
of the universe at that redshift, so that r(m) < tsF is not fulfilled. However, high energy photons 
come mainly from massive stars, which satisfy r(m) < tsp- So at high energies, "Approximation 1" 
results agree well with the full analytical solution. The "Hybrid" approximation however, is more 
accurate through the entire energy range shown in Figure El yet, the results still deviate from the 
full analytical solution. 

Still, the emissivity calculated from the template of Starburst99 is higher than the results 
obtained from the full analytical solution of Eqs. ()AmA4p with fitting formulae. This is mainly 
due to the full stellar evolutionary tracks used by Starburst99, which extend beyond the zero-age 
main sequence (ZAMS) stage. For example, the luminosity of a 7 Mq star with metallicity 1/50 Zq 
at the end of the main sequence is three times as large as the ZAMS luminosity; at the end of its 
evolution, the luminosity is ~ lOx higher compared to ZAMS luminosity. The analogous value for 
a 100 Mq star of the same metallicity, is about 2x the ZAMS luminosity. 
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Fig. 4. — NIRB flux in the 1.25 p, m to 4.5 p,m waveleng th range. The soUd hne is the flux from all 
galaxies from the default model in lHelgason et al.l (|2012l ). while the dashed line is the remaining flux 
after removal of all sources brighter than = 28. The grey regions refer t o the flux range between 
the maximum and minimum evolution scenarios in iHelgason et al.l (|2012l ). As a comparison, we 
plot the flux from all galaxies with z > 5 in our work (crosses) , and the flux after removal of sources 
down to miim = 28 (diamonds). Before any galaxy reni o yah t he flux from galaxies with z > 5 is 
only a few percent of the overall flux in IHelgason et al.l ()2012l ). i.e. the low-z galaxies dominate. 
However, a fter removal of galaxi es with < 28, the flux from the remaining galaxies at all 

redshifts in lHelgason et al.l (|2012l ) is comparable with that from the remaining galaxies with z > 5 
in our work. 
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Fig. 5. — NIRB fluctuations angular power spectrum at different (observer frame) wavelengths, as 
labeled in each panel. The solid (dashed) lines refer to galaxy clustering (shot noise), the thick 
(thin) lines show the case with (without) reni oval of br ight galaxies. We also p lot the observations 
at wavelength 1.6 um (iThompson et al. I l2007al l. 2.4 (jMatsumoto et al.ll201lh . 3.6 fim and 4.5 fim 
([Kashlinskv et al.ll2012l . uncorrected for masking effects) by filled circles with errorbars. 
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Fig. 6.— 



The ratio of 5F/F as a function of I for four different wavelength 
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Fig. 7. — The spectrum of the NIRB fluctuations (sohd line) at / = 2000, the dashed line shows a 
An "'^■'^ law. 
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Fig. 8. — Emissivity of Pop II stars at redshift 10.0. The solid line is the result obtained by using 
the spectrum template from Starburst99 used here, the dotted line is for the "Approximation 1" ; 
dash-dotted line refers to "Approximation 2" . The dash-dotted-dotted-dotted line corresponds to 
the "Hybrid" approximation; finally the dashed line refers to the full analytical solution of Eqs 
()AlllA4p using fitting formulae (see text). 



